Processing seismic data

ABSTRACT

A method of determining an integral of a function over an n-dimensional integration domain from the values of the function at a plurality of discrete data points within the integration domain, comprises partitioning the integration domain into a plurality of k-dimensional simplexes, where k≦n. The function is integrated over each simplex, and the results are summed. The method may be applied in the processing of geophysical data such as, for example, seismic data, where the function is a function of at least the geophysical data. The processed data may then be used to obtain information about the physical properties of the earth&#39;s interior.

FIELD OF THE INVENTION

The present invention relates to a method of processing seismic data orother geophysical data such as electromagnetic data. It is applicable toseismic data processing and migration techniques such as, for example,Kirchhoff migration, beam migration, wave equation migration, multipleattenuation and AVO (amplitude versus offset) analysis.

BACKGROUND OF THE INVENTION

In seismic imaging, as well as in many other areas of seismic dataprocessing, one often needs to compute integrals of the general form:$\begin{matrix}{\int\limits_{D}{A\quad{\mathbb{d}V}}} & (1)\end{matrix}$

In equation (1), A is a function to be integrated, and dV is a volumeelement on D. The integration domain D is a k-dimensional sub-domain ofn-dimensional space (k≦n).

The integration can be over, for example, source and receiver positions,midpoints only, midpoint and offsets or scattering angles etc. Theintegration domain typically is a two-dimensional domain, but it couldbe of higher dimension. The integrand A contains the seismic data, orpart of the seismic data, and may also contain some theoretical terms.Integrals of this type occur in seismic data processing in, for example,migration algorithms (pre stack and post stack migration; time and depthmigration; Kirchhoff, beam and wave equation migration), attenuation ofmultiple reflections in marine data, AVO analysis etc. Where k<n thetopography of the acquisition geometry is taken into account—thishappens when not only the x and y coordinates of the sources andreceivers, but also their z coordinates, are known.

One particular problem in processing seismic data is that the seismicdata, and hence the value of the function A, are generally known only atn_(p) n-dimensional points in the integration domain D. These points atwhich the seismic data are known will hereinafter be referred to as“data points”. Integrals having the form of equation (1) havetraditionally been computed using a ‘binning’ algorithm as described in,for example, “Seismic data processing” by O. Yilmaz in Volume II,Society of Exploration Geophysicists, Tulsa, Okla. p 1019 (2001). Abinning algorithm can be implemented in various ways. For example, onecan evaluate equation (1) by using a summation in which all the datapoints are given equal weights: $\begin{matrix}{{\int_{D}^{\quad}{A{\mathbb{d}V}}}\quad = {{{Vol}(D)}{\sum\limits_{i = 1}^{n_{p}}\quad{A_{i}.}}}} & (2)\end{matrix}$

Vol(D) can be estimated roughly, but is often left out altogether inwhich case equation (2) reduces to: $\begin{matrix}{{\int_{D}^{\quad}{A{\mathbb{d}V}}}\quad = {\sum\limits_{i = 1}^{n_{p}}\quad A_{i}}} & (3)\end{matrix}$

Alternatively one can divide the integration domain D up in a largenumber of ‘bins’ (typically squares or rectangles in the case of atwo-dimensional domain), after which the average value of the function Ain each bin is computed (see e.g. Yilmaz, 2001). The integral is thenevaluated as the sum of all these average values: $\begin{matrix}{{\int_{D}^{\quad}{A{\mathbb{d}V}}}\quad = {\sum\limits_{i}\quad{{\overset{\sim}{A}}_{i}{{{Vol}\left( b_{i} \right)}.}}}} & (4)\end{matrix}$

In equation (4) Ã_(i) denotes the average of all measured values of A inthe i^(th) bin, Vol(b_(i)) is the volume of the i^(th) bin and thesummation is over all the bins.

These prior art binning techniques can work reasonably well if the datapoints are distributed on a regular grid. However, the results aredependent on the binning grid size so that, even if the data points aredistributed on a regular grid, errors can result if an unsatisfactorybinning grid is applied. Furthermore, the data points acquired in aseismic survey are often not distributed over a regular grid—firstly, itis hard in practice to position the seismic sources and receiversexactly on a regular grid and, secondly, the actual optimal survey gridmay be irregular rather than regular. If the data points are irregularlydistributed over the integration domain then the prior art binningtechniques are subject to further errors and become unreliable.

SUMMARY OF THE INVENTION

A first aspect of the present invention provides a method of processinggeophysical data, the method comprising determining an integral of afunction of the geophysical data over an n-dimensional integrationdomain from the values of the function at a plurality of discrete datapoints within the integration domain, the method comprising the stepsof: (a) partitioning the integration domain into a plurality ofk-dimensional simplexes, where k≦n; (b) integrating the function overeach simplex; and (c) summing the results of step (b).

The term “simplex” as used herein denotes the results of anytriangulisation, or its corresponding Voronoi diagrams, of a surface intwo dimensions, a tetrahedron in three dimensions and, in general, ann-dimensional shape having n+1 vertices.

The method of the invention provides more reliable results than theprior art binning techniques. The invention is particularly advantageouswhen applied to a set of irregularly-distributed data points.

A second aspect of the invention provides a method of processinggeophysical data, the method comprising determining an integral of afunction of the geophysical data over an n-dimensional integrationdomain from the values of the function at a plurality of discrete datapoints within the integration domain, wherein the method comprises thesteps of: (a) assigning a weight to each data point; (b) multiplying thevalue of the function at each data point by the weight of the respectivedata point; and (c) summing the results of step (b).

A third aspect of the invention provides a method of processinggeophysical data, the method comprising regularising an array of datapoints located in an n-dimensional domain, each data point representingthe value of a function of the geophysical data at a respective discretelocation within the domain, wherein the method comprises the steps of:(a) partitioning the domain into a plurality of k-dimensional simplexes,where k≦n, with each vertex of a simplex being coincident with arespective one of the data points; (b) determining the co-ordinates of alocation at which it is desired to determine the value of the function,the location being within the domain and not being coincident with oneof the data points; and (c) determining the value of the function at thelocation from the values of the function at the vertices of the simplexthat contains the location.

The invention also provides methods of determining an integral of afunction corresponding to the first and second aspects, and provides amethod of regularising an array of data points corresponding to thethird aspect.

A fourth aspect of the invention provides a method of seismic surveyingcomprising: propagating an acoustic or electromagnetic field through atleast one subsurface layer of the earth; acquiring geophysical data at aplurality of discrete locations; processing the geophysical dataaccording to a method of the first, second or third aspect; anddetermining one or more parameters relating to physical properties ofthe at least one subsurface layer using the processed data.

A fifth aspect of the invention provides an apparatus for processinggeophysical data and adapted to determine an integral of a function ofthe geophysical data over an n-dimensional integration domain from thevalues of the function at a plurality of discrete data points within theintegration domain, the apparatus comprising: (a) means for partitioningthe integration domain into a plurality of k-dimensional simplexes,where k≦n; (b) means for integrating the function over each simplex; and(c) means for summing the results of step (b).

A sixth aspect of the invention provides an apparatus for processinggeophysical data and adapted to determine an integral of a function ofthe geophysical data over an n-dimensional integration domain from thevalues of the function at a plurality of discrete data points within theintegration domain, the apparatus comprising: (a) means for assigning aweight to each data point; (b) means for multiplying the value of thefunction at each data point by the weight of the respective data point;and (c) means for summing the results of step (b).

A seventh aspect of the invention provides an apparatus for processinggeophysical data and adapted to regularise an array of data pointslocated in an n-dimensional domain, each data point representing thevalue of a function of the geophysical data at a respective discretelocation within the domain, the apparatus comprising: (a) means forpartitioning the domain into a plurality of k-dimensional simplexes,where k≦n, with each vertex of a simplex being coincident with arespective one of the data points; (b) means for determining theco-ordinates of a location at which it is desired to determine the valueof the function, the location being within the domain and not beingcoincident with one of the data points; and (c) means for determiningthe value of the function at the location from the values of thefunction at the vertices of the simplex that contains the location.

The apparatus may comprise a programmable data processor.

The invention further provides a storage medium containing a program forthe data processor of such an apparatus.

The invention further provides a storage medium containing a program forcontrolling a programmable data processor to carry out a method of thefirst, second, third or fourth aspect.

BRIEF DESCRIPTION OF THE FIGURES

Preferred embodiments of the present invention will now be described byway of illustrative example with reference to the accompanying figuresin which:

FIG. 1 illustrates tessellation of an integration domain containingirregularly-spaced data points;

FIG. 2(a) illustrates a seismic receiver distributions and a binninggrid coincident with the receiver positions;

FIG. 2(b) illustrates the result of processing seismic data acquired atthe receivers of FIG. 2(a) using the binning grid of FIG. 2(a);

FIG. 3(a) illustrates a seismic receiver distributions and a binninggrid not coincident with the receiver positions;

FIG. 3(b) illustrates the difference between the result of processingseismic data acquired at the receivers of FIG. 3(a) using the binninggrid of FIG. 3(a) and the result of processing the data using a methodof the present invention;

FIG. 4(a) illustrates an irregular seismic receiver distributions and abinning grid not coincident with the receiver positions;

FIG. 4(b) shows the result of processing simulated seismic data usingthe binning grid of FIG. 4(a);

FIG. 4(c) illustrates the result of processing simulated seismic datausing a method of the present invention;

FIG. 5(a) illustrates a random seismic receiver distribution and abinning grid not coincident with the receiver positions;

FIG. 5(b) shows the result of processing simulated seismic data usingthe binning grid of FIG. 5(a);

FIG. 5(c) illustrates the result of processing simulated seismic datausing a method of the present invention;

FIG. 6 illustrates a typical marine seismic surveying arrangement;

FIG. 7 illustrates scattering angles;

FIGS. 8(a) and 8(b) show integration domains in the midpoint/offsetdomain and the scattering angle domain;

FIG. 8(c) shows estimates of the reflectivity determined using a methodof the present invention and three prior art binning techniques;

FIG. 9(a) show an irregular distribution of data points in anintegration domain;

FIG. 9(b) shows estimates of the reflectivity determined by a method ofthe present invention and three prior art binning techniques for thedata points of FIG. 9(a);

FIG. 10(a) show an irregular distribution of data points in anintegration domain;

FIG. 10(b) shows estimates of the reflectivity determined by a method ofthe present invention and three prior art binning techniques for thedata points of FIG. 10(a);

FIG. 11(a) show an irregular distribution of data points in anintegration domain;

FIG. 11(b) shows estimates of the reflectivity determined by a method ofthe present invention and three prior art binning techniques for thedata points of FIG. 11(a);

FIG. 12(a) show an irregular distribution of data points in anintegration domain;

FIG. 12(b) shows estimates of the reflectivity determined by a method ofthe present invention and three prior art binning techniques for thedata points of FIG. 12(a);

FIG. 13 is a block flow diagram of a method of the present invention;

FIG. 14 is a block flow diagram of a regularisation method of thepresent invention; and

FIG. 15 is a block schematic diagram of an apparatus according to theinvention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

In essence, given the positions of the data points the integrationdomain is tessellated into elements T_(i) that constitute simplexes (andso the integration domain would be triangulated in the case of atwo-dimensional integration domain, tessellated into tetrahedrons in thecase of a three-dimensional integration domain, etc). The integration ofthe function A over the integration domain D then reduces to anintegration of the function A over each simplex separately. Thesecontributions then are summed. Thus: $\begin{matrix}{{\int_{D}^{\quad}{A{\mathbb{d}V}}}\quad = {\sum\limits_{i = 1}^{n_{t}}\quad{\int_{T_{i}}^{\quad}{A\quad{\mathbb{d}V_{T_{i}}}}}}} & (5)\end{matrix}$and the summation over i is from i=1 to i=n_(t) where n_(t) is thenumber of simplex elements created in the tessellation.

The integral of the function A over an element can be computed inseveral different ways. The simplest way is to approximate A as aconstant over an element.

Alternatively one can approximate A by a linear function, and such afunction is uniquely determined by its values at the vertices of theelement where the volume element V_(T) _(i) is a k-dimensional simplex(e.g. a triangle (k=2), a tetrahedron (k=3), etc.). It is shown inappendix I below that, for an element T with l (where l=k+1) verticesthe following holds: $\begin{matrix}{{\int_{T}^{\quad}{A{\mathbb{d}V_{T}}}}\quad = {\frac{{Vol}(T)}{k}{\sum\limits_{j = 1}^{l}\quad A_{j}}}} & (6)\end{matrix}$

A_(j) represents the value of the function A at each vertex j of thevolume element V_(T), and the summation is over all vertices (l verticesin total) of the element V_(T). If there is no topography (that is, whenk=n, where n is the dimension of the integration domain) then Vol(T) isthe volume of the element. If there is topography then the expressionfor Vol(T) is slightly more complicated, and is given in appendix II.

Thus the integral of equation (1) has been rewritten as two summations:a summation over j representing summation over the vertices of anelement (the summation is from j=1 to j=l), followed by a summation overall the elements: $\begin{matrix}{{\int_{D}^{\quad}{A{\mathbb{d}V_{D}}}}\quad = {\sum\limits_{i = 1}^{n_{i}}\left( {\frac{{Vol}\left( T_{i} \right)}{k}{\sum\limits_{j = 1}^{l}\quad A_{j}}} \right)}} & (7)\end{matrix}$

It is often more practical to reverse the order of the summation inequation (7) in which case equation (1) may be re-written as:$\begin{matrix}{{\int_{D}^{\quad}{A{\mathbb{d}V_{D}}}}\quad = {\sum\limits_{i = 1}^{n_{p}}{w_{i}A_{i}}}} & (8)\end{matrix}$where w_(i) is called the “weight” of point i and A_(i) is the value ofthe function A for the i^(th) data point. The weight of a data point isdetermined by the areas or volumes of all simplexes that have a vertexat that point.

It is instructive to compare equation (7) with equation (3): thesummation of equation (3) that gives equal weights to all data pointshas been replaced in equation (7) by a summation in which the weightsare determined by the way in which the data points are distributed inthe integration domain. As will be shown in the examples below, theweight of a certain point is relatively small if many other pointssurround it, and it is relatively large if it is an isolated point.

If the data points are distributed on a regular grid then each point hasthe same weight as every other point, and equation (7) becomes identicalto equation (3) or equation (2).

FIG. 13 is a block schematic flow diagram illustrating the inventionapplied to a method of processing geophysical data. Initially at step 1geophysical data such as seismic data or electromagnetic data areacquired. The data are acquired at an array of discrete data points thatare not arranged on a regular grid, so that values for the data areknown only for the data points. Alternatively, the invention may beapplied to pre-existing data, in which case step 1 of acquiring the datais replaced by step 2 of retrieving existing data from storage.

At step 3, the domain in which the data points lie is tessellated. Asexplained above, the domain is tessellated into simplexes, with eachvertex of each simplex being coincident with a respective one of thedata points.

At step 4 the weight of a data point, denoted in FIG. 13 as the i^(th)data point, is determined. Step 4 involves determining the weight ofeach simplex having the i^(th) data point as a vertex.

Steps 3 and 4 are described in more detail with reference to FIG. 1below.

At step 5 the value of the data at the data point is multiplied by theweight determined at step 4. This gives the product w_(i)A_(i) ofequation (8) for that data point.

At step 6 it is determined whether steps 4 and 5 have been carried outfor all data points.

If this step provides a “no” determination, steps 4 and 5 are repeatedfor other data points until a “yes” determination is obtained.

Once a “yes” determination is obtained at step 6, the productsw_(i)A_(i) for each data point are summed at step 7.

The method may be performed in different ways and is not limited to theembodiment of FIG. 13. For example, the weights of every data point maybe determined in a single step rather than in separate steps as in FIG.13. Also, if two or more sets of seismic data are acquired for the samearray of data points it is only necessary to tessellate the domain once;once the domain has been tessellated and the weights of the data pointsdetermined, the weights can be used in the processing of a subsequentset of data acquired at the same array of data points.

FIG. 1 shows a small number of data points that are irregularlydistributed over an integration domain. The data points may represent,for example, midpoints between a seismic receiver and a seismic sourceand so represent points for which seismic data have been acquired in aseismic survey. In a case where evaluation of the integral of equation(1) was desired, the integrand A would be known at these data points butwould not be known elsewhere.

FIG. 1 shows a two-dimensional integration domain and accordingly duringthe tessellation step the domain is triangulated to define a pluralityof triangles as shown in FIG. 1. Each vertex of each triangle is at arespective data point. The triangulation may be carried out using, forexample, a Delaunay triangulation as described by S. W. Sloan in “A fastalgorithm for constructing Delaunay triangulations in a plane”, Adv.Eng. Software, vol. 8, pp 34-55 (1987).

From this triangulation the weights for each point are computed from theareas (in this 2-D case) of the triangles having a vertex at that point(in a higher dimensional case, the weight of a point is determined fromthe volumes of the simplexes having a vertex at that point). The weightof a data point is therefore determined by the distribution ofneighbouring data points—where two data points are said to beneighbouring data points if there is a triangle that has both datapoints as vertices. If, for a particular data point, the neighbouringdata points are distant the triangles with a vertex at that data pointwill have large areas, and so the data point will have a large weight.Conversely, if, for a particular data point, the neighbouring datapoints are close the triangles with a vertex at that data point willhave small areas, and so the data point will have a low weight. Thevalues of some of the weights are shown in FIG. 1. As expected, pointsthat are close to other points have a relatively low weight, whereaspoints that are rather isolated have a relatively high weight.

The integral of equation (1) may then be determined according to theinvention by using equation (5). To determine the integral, the integralof the function A over each triangle is computed from the seismic dataavailable for each vertex of the triangle. Since the vertices of thetriangles are chosen to be coincident with data points, seismic datawill always be available for each vertex, and the integral over eachtriangle may be determined by summing the values of the function A foreach vertex of each element, according to the summation over j inequation (7). The integrals over each triangle are then summed.

Alternatively, once the weights for each data point have been determinedthe integration may be carried out according to the technique ofequation (8), in which the value of the seismic data for each data pointis multiplied by the weight of that point and the weighted values of thedata are then summed.

It should be noted that equations (7) and (8) represent two equivalentmethods of obtaining the integral of equation (5). The same result isobtained by using the weights of the data points in equation (8) as byusing the volume of the simplex elements in equation (7).

Examples of the method of the invention will now be described withreference to determining the earth's reflectivity for seismic energy.

In post stack (or zero-offset) imaging, the reflectivity of the earth isgiven by an integral over zero-offset seismic data. In the time domainthis integral can be written, per D. Miller et al in “A new slant onseismic imaging: Migration and integral geometry”, Geophysics, Vol. 52,pp 943-964. (1987), as: $\begin{matrix}{{R\left( {x,y,z} \right)} = {\int_{D}^{\quad}{{a\left( {m_{1},m_{2},x,y,z} \right)}{u\left( {m_{1},m_{2},{t\quad\quad = \quad{{T_{1}\left( {m_{1},m_{2},x,y,z} \right)} +}}} \right.}}}} & (9) \\{\left. \quad{T_{2}\left( {m_{1},m_{2},x,y,z} \right)} \right){\mathbb{d}\quad m_{1}}{\mathbb{d}\quad m_{2}}} & \quad\end{matrix}$where R is the reflectivity, (m₁,m₂) is the source-receiver midpointlocation, a is the theoretical migration amplitude (containing, amongother quantities, the geometrical spreading, the background velocitymodel and the Beylkin determinant) given by the Generalized RadonTransform in this case, T₁ is the travel time from the midpoint to thescattering point (x,y,z), T₂ is the travel time from the scatteringpoint to the midpoint and u denotes the seismic data. Equation (9) is inthe form of equation (1), with the product au of the theoreticalmigration amplitude and the seismic data representing the function A ofequation (1). Equation (9) therefore can be evaluated using a prior artbinning technique (equations (2), (3) or (4)) or using a method of theinvention.

FIGS. 2(a) to 5(c) illustrate the improved results obtained by a methodof the present invention.

FIG. 2(a) illustrates an array of receiver positions, with each receiverposition denoted by an “×”. Synthetic seismic data were generated foreach receiver location shown in FIG. 2(a) using a simple 3-layer modelfor the earth's interior. That is, the earth's interior was representedby three layers of different seismic properties, with the upper twolayers having a finite thickness and with the third layer extending toinfinite depth.

The synthetic seismic data simulated for each receiver location of FIG.2(a) were then processed to obtain information about the earth'sinterior. The data were processed in two ways:—firstly, using the priorart binning technique of equation (4) and, secondly, using a method ofthe invention according to equation (8) above. If the processing isaccurate, the original three-layer earth model should be recovered.

FIG. 2(a) also shows the binning grid applied, in solid lines. It willbe seen that each vertex of the binning grid 1 is coincident with areceiver location, and that each receiver location is coincident with avertex of the binning grid. The binning grid is coincident with thereceiver grid.

FIG. 2(b) illustrates the results of processing the simulated seismicdata. It will be seen that the 3-layer model structure of the earth isrecovered well when the data are processed—three layers 2, 3 and 4 canclearly be discerned in FIG. 2(b).

The receiver locations in FIG. 2(a) are arranged in a regular grid, andthe size of each “bin” of the binning grid is equal to the distancebetween adjacent receiver locations. Accordingly, the weights that areallocated to each data point are the same, and the results of the methodof the invention are identical to the results of the binning technique.The results of FIG. 2(b) are therefore obtained either using the priorart binning technique of equation (4) or using the method of theinvention.

FIGS. 3(a) and 3(b) illustrate the results of using a different binninggrid. FIG. 3(a) again shows the receiver locations, with each receiverlocation being denoted by an “×”. The receiver locations of FIG. 3(a)are the same as the receiver locations of FIG. 2(a), and the seismicdata simulated for the receiver locations of FIG. 3(a) using the simplethree-layer model of the earth are therefore identical to the seismicdata simulated for the receiver locations of FIG. 2(a). However, thebinning grid 5 applied in FIG. 3(a) is not the same as the binning grid1 applied in FIG. 2(a). The binning grid 5 applied in FIG. 3(a) haslarger bins, and the vertices of the binning grid are, in general, notcoincident with receiver locations. The majority of receivers aretherefore not positioned at a vertex of the binning grid. Each bincontains more than one receiver.

FIG. 3(b) shows the difference between, on the one hand, the result ofprocessing the seismic data simulated for the receiver locations of FIG.3(a) according to the prior art technique of equation (4) using thebinning grid 5 of FIG. 3(a) and, of the other hand, the result ofprocessing the simulated seismic data using the method of the inventionaccording to equation (8). As can be seen from FIG. 3(b) there aresignificant differences between the results acquired by the twoprocessing techniques. In FIG. 3(b) the difference is shown as afunction of depth within the earth (vertical scale) and horizontallocation (horizontal scale). Dark areas in FIG. 2(b) represent thegreatest positive difference, where the result obtained by the inventionminus the result obtained by a binning technique is positive, and lightareas represent the greatest negative difference.

The results of processing the simulated seismic data according toequation (8) do not depend on the binning grid, so that equation (8)will give the same results when applied to the seismic data simulatedfor FIG. 2(a) as when applied to the seismic data simulated for FIG.3(b). The simulated seismic data of FIG. 2(a) are the same as thesimulated seismic data of FIG. 3(a), as they were simulated for the samereceiver array using the same earth model. The difference between thetwo results shown in FIG. 3(b) is therefore entirely due to a change inthe results obtained using the prior art binning technique.

Furthermore, the difference between the two sets of results in FIG. 3(b)is equivalent to the difference between the results given by the binninggrid 5 of FIG. 3(a) and the results given by the binning grid 1 of FIG.2(a). This indicates that the accuracy of the prior art binningtechnique is heavily dependent on the choice of the binning grid used.In fact, the error in an integral performed by the prior art binningtechnique is related to the number of data points in each bin cell(otherwise known as the “fold”).

FIG. 4(a) illustrates a further array of receiver locations, with eachreceiver location being indicated by an “×”. The receivers are arrangedgenerally in a grid, although the receiver arrangement is not exactlyregular and some receivers are slightly displaced from the position theywould occupy if the grid were entirely regular. Seismic data weresimulated for each receiver location of FIG. 4(a), using the same simple3-layer model for the earth's interior as used in the two previousexamples.

The simulated seismic data were then processed according to the priorart binning technique of equation (4) and according to the method ofequation (8) of the invention. The processing according to equation (4)was carried out using the binning grid 5 shown in FIG. 4(a), which isthe same as the grid 5 used in FIG. 3(a).

FIG. 4(b) shows the difference between the results obtained using thebinning algorithm in this example and the results of FIG. 2(b). FIG.4(c) shows the difference between the results of processing the seismicdata of this example using equation (8) and the results of FIG. 2(b). Asin FIG. 3(b), the difference is plotted against depth in the earth andhorizontal location, and the magnitude of the difference is representedby shading.

As noted above, the results of FIG. 2(b) are essentially correct, sothat the results of FIG. 4(b) represent the error in processing thesimulated seismic data for the receiver locations of FIG. 4(a) using thebinning technique, and the results of FIG. 4(c) represent the error inprocessing the seismic data for the receiver location shown in FIG. 4(a)using the method of the present invention. It will be seen from FIG.4(c) that the error in the results obtained by the method of the presentinvention is remarkably small. However, the results of FIG. 4(b) showthat the prior art binning technique provides much less accurate resultsfor the slightly irregular receiver array of FIG. 4(a).

FIG. 5(a) shows a further array of receivers. The receiver positions,which are again denoted by an “×”, were determined randomly, so that thearrangement of receivers is completely irregular.

Seismic data were simulated for the receivers shown in FIG. 5(a) usingthe simple 3-layer earth model used in previous examples. The simulatedseismic data were then processed according to the prior art binningtechnique of equation (4) using the binning grid 5 of the two previousexamples, and were also processed according to equation (8). FIG. 5(b)shows the error in the results obtained using the prior art binningtechnique, and FIG. 5(c) shows the errors in the results obtained usingthe method of the present invention. As in FIGS. 4(b) and 4(c), theerrors in the results of FIGS. 5(b) and 5(c) are defined as being thedifference between the results and the results of FIG. 2(b). The resultsare plotted in the same manner as the results of FIGS. 3(b), 4(b) and4(c).

It will again be seen that the errors in the method of the presentinvention (FIG. 5(c)) are significantly lower than the errors in theprior art binning technique. This illustrates that the method of thepresent invention is much more reliable and can be applied far morebroadly than the prior art binning technique. The prior art binningtechnique can provide good results if the data points are spacedregularly, and if the binning grid is chosen carefully. However if thebinning grid is not chosen carefully, or if the data points are notregularly arranged, the prior art binning technique does not give goodresults. The present invention provides much better results in the caseof an irregular array of data points; furthermore, the present inventiondoes not require a binning grid, so the possibility of introducingerrors through use of an inappropriate binning grid are eliminated.

The above example has been described with reference to the receiverpositions, and this is appropriate for zero-offset imaging. Theinvention may alternatively be effected using, for example,source-receiver midpoint locations, and this would be appropriate forpost-stack migration.

FIGS. 6 to 12 illustrate a further example of the present invention.This example relates to seismic data acquired using a seismic cable asshown in FIG. 6. In this seismic surveying arrangement, a plurality ofseismic receivers (not shown) are mounted on a receiver cable 6. Thereceiver cable 6 has a typical length of a few kilometres, for example 3km. The receivers are uniformly spaced along the cable 6, with thespacing between two neighbouring receivers being, for example, 12.5 m or25 m. A seismic source 7 is also mounted on the seismic cable 6, and ispositioned substantially midway along the length of the cable. Thecable, with the source, is moved along the earth's surface, and thesource 7 is repeatedly “fired” (that is, is repeatedly actuated to emitseismic energy into the earth). The time interval between successivefirings of the source is controlled so that the source movesapproximately 50 m between one firing position and the next.

The model of the earth's interior assumes that there is a singlereflector 8 of seismic energy, which is substantially horizontal.

Seismic data were again simulated for the seismic surveying arrangementof FIG. 6, for the simple earth model having the single flat reflector 8of FIG. 6.

The data were processed to obtain the reflectivity R by applying atwo-dimensional pre-stack Kirchhoff integral according to:R(x,z)=∫a(θ,φ,x,z)u(θ,φ,t=T ₁(θ,φ,x,z)+T ₂(θ,φ,x,z))dθdφ  (10)

In equation (10), as in equation (9), the function a contains atheoretical amplitude. The variable u represents the seismic data, whichare known only at the data points. θ,φ are scattering angles, and aredefined in FIG. 7. The quantities p_(r), p_(s) defined in FIG. 7 are theslowness of the ray from the receiver to the scattering point (x,z), andthe slowness of the ray from the source to the scattering point. (Theslowness is the gradient of the travel time function). The quantityp_(r)+p_(s) defined in FIG. 7 is the sum of p_(r) and p_(s). Equation(10) is a two-dimensional equation, but it may be extended tothree-dimensions as proposed by S Brandsberg-Dahl et al in “Focussing inDip and AVA Compensation on Scattering Angle/Azimuth Common ImageGather”, Geophysics, Vol. 68, pp 232-254 (2003).

The use of scattering angles as integration variables in equation (10),as opposed to using midpoint and offset as the integration variables asin equation (9) has a number of advantages. There is no need to computeand store the Beylkin determinant, and multi-pathing can beautomatically included. Possible disadvantages are that the integrationdomain becomes more irregular, and that the triangulation changes fromone scattering point to the next, thereby increasing the computationaleffort required to determine the integration. However, the irregularityof the integration domain can be overcome by the use of equation (8),involving the weights of the data points. The problem of thetriangulation changing from one scattering point to the next can beovercome by performing triangulation in the midpoint-offset domain—thereis a one-to-one relationship between the midpoint-offset domain and thescattering angle domain with a strictly positive Jacobian (assumingthere is no multi-pathing) so that triangulation need be performed onlyonce. Thus, the method of equation (7) may also be used, by performingtriangulation in the midpoint-offset domain and then performing theremainder of the integration in the scattering angle domain—if theintegration is done in this way, there is no additional computationaleffort required. If there is multi-pathing, it may still be possible toperform a triangulation in the midpoint-offset domain, provided that theresultant triangulation is modified to account for multi-pathing, andthis allows equation (7) to be used with negligible additionalcomputational costs.

The synthetic seismic data simulated for the acquisition arrangement ofFIG. 6 were then processed according to a prior art binning techniqueand according to a method of the present invention. The expected resultof processing the simulated seismic data is, of course, that the flatreflector earth structure of FIG. 6 should be recovered.

FIG. 8(c) shows results obtained using an acquisition geometry that isregular in the midpoint-offset domain. The source spacing was exactly12.5 m, and the receiver spacing was exactly 25 m. This leads to theregular array of data points for the midpoint offset domain as shown inFIG. 8(a). FIG. 8(b) shows the corresponding array of data points in thescattering angle domain, and it will be seen that the array of datapoints is not completely regular in the scattering angle domain.

FIG. 8(c) shows the reflectivity obtained by processing the seismic datasimulated for the data points of FIG. 8(a). The solid line in FIG. 8(c)shows the reflectivity obtained using the method of equation (8).

The seismic data were also processed using the prior art binning methodof equation (4). This method was performed using three different binninggrids, of three degrees (results for this binning grid are shown in FIG.8(c) as a dotted line), six degrees (the results for this binning gridare shown in FIG. 8(c) as a dot-dash line), and nine degrees (theresults for this binning grid are shown in FIG. 8(c) as a dashed line).

It will be seen that the reflectivity for the method of the invention,and the reflectivities obtained by the prior art binning technique forall three binning grids are very close to one another.

FIG. 9(a) shows an array of data points obtained by randomly perturbingthe source position and the cable position each by 50 m. The data pointsare shown in the scattering angle domain in FIG. 9(a), and it will beseen that the array is slightly irregular.

FIG. 9(b) shows the results of processing seismic data simulated for thedata points of FIG. 9(a) using the simple earth model of FIG. 6. FIG.9(b) shows the reflectivity obtained by use of the method of equation(8) and by use of the prior art method of equation (4) for the threedifferent binning grids used in the example of FIGS. 8(a) to 8(c). Thesolid line in FIG. 9(b) represents the results of equation (8), and thedotted line, the dot-dashed line, and the dashed line in FIG. 9(b)represent the results of the three prior art binning techniques. Thedotted line, dot-dashed line and dashed lines in FIG. 9(b) denote thesame binning grids as the dotted line, dot-dashed line, and dashed lineof FIG. 8(c), and this is also the case for the results in FIGS. 10(b),11(b) and 12(b) below.

FIG. 9(b) clearly shows that the irregularities in the array of datapoints has a marked effect on the results obtained by the prior artbinning technique. The prior art binning technique does not obtain aflat reflectivity, for any of the three bin sizes used. However, thereflectivity obtained by equation (8) is still substantially flat, andis much flatter than the results given by the prior art binningtechniques. It is clear that the method of equation (8) provides asignificantly more accurate result than the prior art binningtechniques.

FIG. 10(a) shows an array of data points in the scattering angle domainthat corresponds generally to the array of FIG. 9 but with 5% of thereceivers removed. The results of processing seismic data simulated forthe data points of FIG. 10(a) using the simple earth model of FIG. 6 areshown in FIG. 10(b). FIG. 10(b) shows the reflectivity obtained usingthe method of equation (8) and using the binning technique using thethree binning grids of the previous two examples. It will again be seenthat the reflectivity results given by the binning techniques havedeteriorated slightly compared with the reflectivity results obtained bythe binning technique in FIG. 9(b). However, the method of equation (8)still yields a reflectivity that is substantially flat, and that is farmore accurate than the reflectivity given by any of the binningtechniques.

FIG. 11 shows an array of data points in the scattering angle domain,for a receiver spacing of 25 m. The source spacing was still 12.5 m. Thecorresponding reflectivity results obtained from seismic data simulatedfor these data points (using the simple earth model of FIG. 6) are shownin FIG. 11(b). FIG. 11(b) shows the reflectivity obtained using themethod of equation (8) and using the binning technique using the threebinning grids of the previous two examples. If the results are comparedwith those of FIG. 8(b), it will be seen that the method of equation (8)still gives a substantially flat reflectivity. The reflectivity given bythe binning algorithms is, however, noticeably less flat than thereflectivity given in FIG. 8(b). This shows that increasing the receiverspacing, so that the array of data points is less dense, increases theerror in results obtained using a conventional binning technique. Themethod of the present invention is, however, unaffected.

FIG. 12(a) shows a further array of data points in the scattering angledomain. This array was again obtained using a nominal receiver spacingof 25 m, although the source and receiver positions have been perturbedand 5% of the receivers have been removed. This yields an irregulararray of data points as clearly shown in FIG. 12(a).

Seismic data were simulated for the data points of FIG. 12(a) using thesimple one-layer model of FIG. 6. FIG. 12(b) shows the reflectivityobtained by processing the simulated seismic data, and shows thereflectivity obtained using the method of equation (8) and using theprior art binning technique for the three binning grids of the previousfour examples. It will be seen that the results obtained by the binningtechnique are dependent heavily on the binning grid used—that is, theresults obtained using the prior art binning techniques are veryunstable. Furthermore, none of the prior art binning techniques produceda reflectivity that is substantially flat. As in the previous figures,the method of the present invention again gives the best results for thereflectivity.

The results of FIGS. 8(a) to 12(b) show clearly that the method of thepresent invention provides far more reliable results than the prior artbinning techniques. The method of the invention provides good resultsregardless of whether the array of data points is regular or irregular.In contrast, a prior art binning technique will provide good resultsonly for a regular array of data points, and only if the binning grid ischosen correctly. In the case of an irregular array of data points, theresults given by a prior art binning technique depend heavily on thebinning grid used.

Integrals of the form of equation (1) are also used in electromagneticimaging of the earth. Electromagnetic imaging has important applicationsin determining the earth's structure in the vicinity of the borehole(see, for example, Y. Chen and M. Oristaglio in “A modelling study ofborehole radar for oilfield applications”, Geophysics, 67, 1486-1494(2002).

In electromagnetic imaging the acquired data are the values of theelectric and magnetic fields and the goal is to determine the electricalconductivity and permittivity of the earth. For example the electricconductivity is given (see equation (30) of T. Wang and M. Oristaglio in“GPR imaging using the generalised Radon Transform”, Geophysics Vol. 65,1553-1559 (2000)) as: $\begin{matrix}{{\delta\quad{\sigma(x)}} = {\int{{f\left( {E,B} \right)}{\mathbb{d}\alpha}{\mathbb{d}\beta}}}} & (11)\end{matrix}$where δσ(x) is the electric conductivity at a point x, f is a knownfunction, E and B are the electric and magnetic fields respectively, andα and β are angles. In practice, the values of E and B will be knownonly at a finite number of points. If these points are distributed on anirregular grid, then equation (11) may be determined using a method ofthe invention. For example equation (11) may be re-written along thelines of equation (7) or equation (8), and the same algorithms as in theprevious examples may be applied.

The integration methods described so far use the exact source andreceiver locations.

They also assume an appropriate background model for, for example, thepropagation velocity of seismic energy through the earth. In practicethere is an uncertainty in all these parameters—for example, thepositions of the seismic sources and seismic receivers will not beprecisely know. If it is assumed that these uncertainties can bequantified using a statistical distribution, such as for example aGaussian distribution, then it is also possible to estimate theuncertainty in the integral itself using the methods of the inventiondescribed above.

Thus, for example, equation (10) can be rewritten as: $\begin{matrix}{< {R\left( {x,z} \right)}>={\int_{C}^{\quad}{< {a\left( {x,z,\theta,\varphi} \right)} > {u\left( {\theta,\varphi,t} \right.}}}} & (12) \\{\left. \quad{= {< {{T_{1}\left( {\theta,\varphi,x,z} \right)} + {T_{2}\left( {\theta,\varphi,x,z} \right)}} >}} \right){\mathbb{d}\theta}{\mathbb{d}\varphi}} & \quad\end{matrix}$where <q> denotes the value of the parameter q with a certainstatistical distribution.

Note that the amplitude a and travel times T₁ and T₂ implicitly dependon the source and receiver positions and the velocity model.

The value of <R> can be determined by finding a range of values for thesource and receiver positions and for the velocity model, sampledaccording to their distributions.

This gives a range of values for the amplitude a and travel times T₁ andT₂, which in turn, after computation of the integral of equation (12),gives a range of values for R. A distribution for R can be obtained fromthis range of value for R, and this distribution provides an estimate ofthe uncertainty in R.

The invention has been described above with reference to processinggeophysical data, but the invention has further applications. Oneapplication is in regularising a set of data points.

As explained above, the data points obtained in an actual seismic surveyare often not arranged on a regular grid. In this case, in aconventional processing method the binning step is frequently precededby a regularisation step. In the regularisation step values for theintegrand A for an array of estimated data points that are on a regulargrid are obtained from the values known for the irregularly-distributeddata points, by interpolation.

The present invention may be used to perform the regularisation step,for example in the processing of geophysical data, which may beperformed very efficiently once the integration domain has beentessellated. This is illustrated in FIG. 14, which is a block schematicdiagram of a method of using the techniques of the invention toperforming a regularisation step. Initially at step 1 geophysical datasuch as seismic data or electromagnetic data are acquired for an arrayof data points that are not arranged on a regular grid. Alternatively,the invention may be applied to pre-existing data, in which case step 1of acquiring the data is replaced by step 2 of retrieving existing datafrom storage.

At step 3, one or more points at which it is desired to estimate valuesfor the seismic data are defined. Step 3 may consist of, for example,defining an array of points that lie on a regular grid.

At step 4, the domain in which the data points lie is tessellated. Theorder in which steps 3 and 4 are performed is not important, and step 4may alternatively be performed before, or at the same time as, step 3.

At step 5 one of the points defined at step 3 is selected, and at step 6it is determined whether this point is co-incident with a data point. Ifthe selected point is not coincident with a data point, the simplexwithin which the selected data point lies is determined at step 7. Atstep 8 the value of the data at the selected point is estimated fromvalues of the data at the vertices of the simplex identified in step 8,for example using linear interpolation.

If step 6 determines that the selected point is coincident with one ofthe data points, the value of the data for that data point is taken tobe the value of the data for the selected point without any furthercalculation being required (this step is omitted from FIG. 14).

FIG. 15 is a schematic block diagram of a programmable apparatus 10according to the present invention. The apparatus comprises aprogrammable data process 11 with a programme memory 12, for instance inthe form of a read-only memory (ROM), storing a programme forcontrolling the data processor 11 to perform any of the processingmethods described above. The apparatus further comprises non-volatileread/write memory 13 for storing, for example, any data which must beretained in the absence of power supply. A “working” or scratch padmemory for the data processor is provided by a random access memory(RAM) 14. An input interface 15 is provided, for instance for receivingcommands and data. An output interface 16 is provided, for instance fordisplaying information relating to the progress and result of themethod. Data for processing may be supplied via the input interface 14,or may alternatively be retrieved from a machine-readable data store 17.

The programme for operating the system and for performing any of themethods described hereinbefore is stored in the programme memory 12,which may be embodied as a semi-conductor memory, for instance of thewell-known ROM type. However, the programme may be stored in any othersuitable storage medium, such as magnetic data carrier 12 a, such as a“floppy disk” or CD-ROM 12 b.

Appendix I

In this appendix we provide proof equation (6) for two and threedimensions. The generalization to higher dimensions follows readily.

Let D be an area in the two dimensional plane. If the area of D is equalto A then by definition A = ∫_(D)    𝕕x𝕕y.

We compute the integral of a linear function over D: $\begin{matrix}{{{\int{\left( {c_{0} + {c_{1}x} + {c_{2}y}} \right){\mathbb{d}x}{\mathbb{d}y}}} = {{{c_{0}{\int_{D}^{\quad}\quad{{\mathbb{d}x}{\mathbb{d}y}}}} + {c_{1}{\int_{D}^{\quad}\quad{x{\mathbb{d}x}{\mathbb{d}y}}}} + {c_{2}{\int_{D}^{\quad}\quad{y{\mathbb{d}x}{\mathbb{d}y}}}}} = {A\left( {c_{0} + {c_{1}\frac{\int_{D}^{\quad}\quad{x{\mathbb{d}x}{\mathbb{d}y}}}{\int_{D}^{\quad}\quad{{\mathbb{d}x}{\mathbb{d}y}}}} + {c_{2}\frac{\int_{D}^{\quad}\quad{y{\mathbb{d}x}{\mathbb{d}y}}}{\int_{D}^{\quad}\quad{{\mathbb{d}x}{\mathbb{d}y}}}}} \right)}}}{{Now}\left( {\frac{\int_{D}^{\quad}\quad{x{\mathbb{d}x}{\mathbb{d}y}}}{\int_{D}^{\quad}\quad{{\mathbb{d}x}{\mathbb{d}y}}},\frac{\int_{D}^{\quad}\quad{y{\mathbb{d}x}{\mathbb{d}y}}}{\int_{D}^{\quad}\quad{{\mathbb{d}x}{\mathbb{d}y}}}} \right)}} & ({A1})\end{matrix}$is the centre of gravity of D. Therefore if D is a triangle with corners(x₁,y₁),(x₂,y₂),(x₃,y₃) then its centre of gravity is${\frac{1}{3}\left( {{x_{1} + x_{2} + x_{3}},{y_{1} + y_{2} + y_{3}}} \right)},$so that: $\begin{matrix}{{\int_{D}^{\quad}{\left( {c_{0} + {c_{1}x} + {c_{2}y}} \right)\quad{\mathbb{d}x}{\mathbb{d}y}}} = {A\left( {c_{0} + {c_{1}\frac{x_{1} + x_{2} + x_{3}}{3}} + {c_{2}\frac{y_{1} + y_{2} + y_{3}}{3}}} \right)}} & ({A2})\end{matrix}$

Or, in words: the integral of a linear function over a triangle is equalto the average of the value of the function at the three vertices of thetriangle.

This result can easily be generalized to higher dimensions. For example,in 3 dimensions, when D is a tetrahedron with volume V and coordinates(x_(i),y_(i),z_(i)) (i=1,2,3,4) we have: $\begin{matrix}{{\int_{D}^{\quad}\quad{\left( {c_{0} + {c_{1}x} + {c_{2}y} + {c_{3}z}} \right){\mathbb{d}x}{\mathbb{d}y}{\mathbb{d}z}}} = {\frac{1}{4}{V\left( {a_{1} + a_{2} + a_{3} + a_{4}} \right)}}} & ({A3})\end{matrix}$where the a_(i) are the values of the linear function at the vertices ofthe tetrahedron.Appendix II—The Volume of a k-dimensional Tetrahedron in n-dimensionalSpace

Consider a k-dimensional tetrahedron T defined by the k+1 n-dimensionalpoints 0,a₁, . . . , a_(k). We will show that the volume of thistetrahedron is given by: $\begin{matrix}{{{Vol}(T)} = {\frac{1}{k!}{\begin{matrix}\left( {a_{1},a_{1}} \right) & \left( {a_{1},a_{2}} \right) & \cdots & \left( {a_{1},a_{k}} \right) \\\left( {a_{2},a_{1}} \right) & ⋰ & \quad & \vdots \\\vdots & \quad & \quad & \quad \\\left( {a_{k},a_{1}} \right) & \cdots & \quad & {\quad\left( {a_{k},a_{k}} \right)}\end{matrix}}^{1/2}}} & ({A4})\end{matrix}$where (.,.) denotes the n-dimensional inner product.

Equation (A4) can be interpreted as:Volume=Base×Height

We prove this using the induction principle. Consider a k-dimensionalparallelepiped in n dimensions spanned by the vectors a₁, . . . , a_(k).The vectors a₁, . . . ,a_(k−1) span the ‘base’ of this parallelepiped.We write:a _(k) =a _(b) +a _(b)where a_(b)=λ₁a₁+ . . . λ_(k−1)a_(k−1)and a_(b) is perpendicular to all a₁, . . . , a_(k−1). |a_(h)| is theheight of the parallelepiped. We now have ${\begin{matrix}\left( {a_{1},a_{1}} \right) & \left( {a_{1},a_{2}} \right) & \cdots & \left( {a_{1},a_{k}} \right) \\\left( {a_{2},a_{1}} \right) & ⋰ & \quad & \vdots \\\vdots & \quad & \quad & \quad \\\left( {a_{k},a_{1}} \right) & \cdots & \quad & {\quad\left( {a_{k},a_{k}} \right)}\end{matrix}} = {{\begin{matrix}\left( {a_{1},a_{1}} \right) & \left( {a_{1},a_{2}} \right) & \cdots & \left( {a_{1},a_{k}} \right) \\\left( {a_{2},a_{1}} \right) & ⋰ & \quad & \vdots \\\vdots & \quad & \quad & \quad \\\left( {a_{k},a_{1}} \right) & \cdots & \quad & {\quad{\left( {a_{b},a_{k}} \right) + \left( {a_{k},a_{k}} \right)}}\end{matrix}} = {\left( {a_{h},a_{h}} \right)\quad{{\begin{matrix}\left( {a_{1},a_{1}} \right) & \left( {a_{1},a_{2}} \right) & \cdots & \left( {a_{1},a_{k}} \right) \\\left( {a_{2},a_{1}} \right) & ⋰ & \quad & \vdots \\\vdots & \quad & {\quad ⋰} & \quad \\\left( {a_{k},a_{1}} \right) & \cdots & \quad & {\quad\left( {a_{k - 1},a_{k}} \right)}\end{matrix}}.}}}$

Upon taking the square root we find equation (A4). Since theparallelepiped consists of k! equal tetrahedrons the volume of thetetrahedron follows.

1. A method of processing geophysical data, the method comprisingdetermining an integral of a function of the geophysical data over ann-dimensional integration domain from the values of the function at aplurality of discrete data points within the integration domain, whereinthe method comprises the steps of: (a) partitioning the integrationdomain into a plurality of k-dimensional simplexes, where k≦n; (b)integrating the function over each simplex; and (c) summing the resultsof step (b).
 2. A method of determining an integral of a function overan n-dimensional integration domain from the values of the function at aplurality of discrete data points within the integration domain, themethod comprising the steps of: (a) partitioning the integration domaininto a plurality of k-dimensional simplexes, where k≦n; (b) integratingthe function over each simplex; and (c) summing the results of step (b).3. A method as claimed in claim 1 or 2 wherein each vertex of a simplexis coincident with a respective one of the data points.
 4. A method asclaimed in claim 1, 2 or 3 and comprising, for a selected simplex,approximating the function as a constant over the simplex.
 5. A methodas claimed in claim 1, 2 or 3 and comprising, for a selected simplex,approximating the function as a linear function over the simplex.
 6. Amethod as claimed in claim 4 or 5 wherein step (b) comprises, for aselected simplex, determining the average of the values of the functionat each vertex of the simplex.
 7. A method of processing geophysicaldata, the method comprising determining an integral of a function of thegeophysical data over an n-dimensional integration domain from thevalues of the function at a plurality of discrete data points within theintegration domain, wherein the method comprises the steps of: (a)assigning a weight to each data point; (b) multiplying the value of thefunction at each data point by the weight of the respective data point;and (c) summing the results of step (b).
 8. A method of determining anintegral of a function over an n-dimensional integration domain from thevalues of the function at a plurality of discrete data points within theintegration domain, the method comprising the steps of: (a) assigning aweight to each data point; (b) multiplying the value of the function ateach data point by the weight of the respective data point; and (c)summing the results of step (b).
 9. A method as claimed in claim 7 or 8wherein step (a) comprises assigning weights to a data point on thebasis of the distribution of neighbouring data points over theintegration domain.
 10. A method as claimed in any of claims 7 to 9wherein the weight assigned to a selected data point is inverselyrelated to the number of other data points in the vicinity of theselected data point.
 11. A method as claimed in any of claims 7 to 10wherein step (a) comprises partitioning the integration domain into aplurality of k-dimensional simplexes, where k≦n, with each vertex of asimplex being coincident with a respective one of the data points, andwherein the weight assigned to a selected data point is proportional tothe volume of the simplex associated with the selected data point.
 12. Amethod of processing geophysical data, the method comprisingregularising an array of data points located in an n-dimensional domain,each data point representing the value of a function of the geophysicaldata at a respective discrete location within the domain, wherein themethod comprises the steps of: (a) partitioning the domain into aplurality of k-dimensional simplexes, where k≦n, with each vertex of asimplex being coincident with a respective one of the data points; (b)determining the co-ordinates of a location at which it is desired todetermine the value of the function, the location being within thedomain and not being coincident with one of the data points; and (c)determining the value of the function at the location from the values ofthe function at the vertices of the simplex that contains the location.13. A method of regularising an array of data points located in ann-dimensional domain, each data point representing the value of afunction at a respective discrete location within the domain, the methodcomprising the steps of: (a) partitioning the domain into a plurality ofk-dimensional simplexes, where k≦n, with each vertex of a simplex beingcoincident with a respective one of the data points; (b) determining theco-ordinates of a location at which it is desired to determine the valueof the function, the location being within the domain and not beingcoincident with one of the data points; and (c) determining the value ofthe function at the location from the values of the function at thevertices of the simplex that contains the location.
 14. A method asclaimed in claim 12 or 13 wherein step (c) comprises determining thevalue of the function at the location by linear interpolation from thevalues of the function at the vertices of the simplex that contains thelocation.
 15. A method as claimed in any preceding claim wherein k=2whereby each simplex is a triangle.
 16. A method as claimed in any ofclaims 1 to 14 wherein k=3 whereby each simplex is a tetrahedron.
 17. Amethod as claimed in any of claims 1, 7 or 12, or in any one of claims 3to 6, 9 to 11, and 14 to 16 when dependent directly or indirectly fromclaim 1, 7 or 12, wherein the function is a function of the geophysicaldata and of at least one further parameter.
 18. A method as claimed inany of claims 1, 7 or 12, or in any one of claims 3 to 6, 9 to 11, and14 to 16 when dependent directly or indirectly from claim 1, 7 or 12,and comprising the further step of determining one or more parametersrelating to physical properties of the earth's interior from theprocessed geophysical data.
 19. A method as claimed in any of claims 2,8 or 13, or in any one of claims 3 to 6, 9 to 11, and 14 to 16 whendependent directly or indirectly from claim 2, 8 or 13, wherein thefunction is a function of geophysical data.
 20. A method of seismicsurveying comprising: propagating an acoustic or electromagnetic fieldthrough at least one subsurface layer of the earth; acquiringgeophysical data at a plurality of discrete locations; processing thegeophysical data according to a method as defined in any one of claims 1to 17; and determining one or more parameters relating to physicalproperties of the at least one subsurface layer using the processeddata.
 21. An apparatus for processing geophysical data and adapted todetermine an integral of a function of the geophysical data over ann-dimensional integration domain from the values of the function at aplurality of discrete data points within the integration domain, theapparatus comprising: (a) means for partitioning the integration domaininto a plurality of k-dimensional simplexes, where k≦n; (b) means forintegrating the function over each simplex; and (c) means for summingthe results of step (b).
 22. An apparatus for processing geophysicaldata and adapted to determine an integral of a function of thegeophysical data over an n-dimensional integration domain from thevalues of the function at a plurality of discrete data points within theintegration domain, the apparatus comprising: (a) means for assigning aweight to each data point; (b) means for multiplying the value of thefunction at each data point by the weight of the respective data point;and (c) means for summing the results of step (b).
 23. An apparatus forprocessing geophysical data and adapted to regularise an array of datapoints located in an n-dimensional domain, each data point representingthe value of a function of the geophysical data at a respective discretelocation within the domain, the apparatus comprising: (a) means forpartitioning the domain into a plurality of k-dimensional simplexes,where k≦n, with each vertex of a simplex being coincident with arespective one of the data points; (b) means for determining theco-ordinates of a location at which it is desired to determine the valueof the function, the location being within the domain and not beingcoincident with one of the data points; and (c) means for determiningthe value of the function at the location from the values of thefunction at the vertices of the simplex that contains the location. 24.An apparatus as claimed in any of claims 21, 22 and 23 and comprising aprogrammable data processor.
 25. A storage medium containing a programfor the data processor of an apparatus as defined in claim
 24. 26. Astorage medium containing a program for controlling a programmable dataprocessor to carry out a method as defined in any of claims 1 to 20.